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MODELING OF ZERO GRAVITY VENTING 


The University of Michigan 
Ann Arbor. MI 


ABSTRACT 


An experimental investigation of the venting of cylindrical containers partially filled 
witn initially saturated liquids was previously conducted under zero-gravity conditions 
at the NASA Lewis Research Center 5-second zero gravity facility, and compared with 
an analytical model which determined the effect of interfacial mass transfer on the 
uliege pressure response during venting. A new model is proposed here to improve the 
estimation of the interfacial mass transfer. Duhammel’s superposition integral is 
incoroorated in this analysis to approximate the transient temperature response of the 
interface, treating the liquid as a semi-infinite solid with conduction heat transfer. 
The results show that this approach to estimating interfacial mass transfer gives 
improvec response when compared to previous models. However, the present model still 
predicts e pressure decrease greater than those in the experiments reported. 


INTRODUCTION 


Tne use of high-energy liquid propellants in the space program has led to a need for 
Information concerning tne thermodynamic bena'ior of cryogenic fluids in tanks which 
are vented or Depressurized to space. Low vapor vent rates are used es a method of 
tank pressure control. The task of venting in low gravity has been successfully 
accomplished curing a number of past missions with venting systems that rely 
exclusively on auxiliary thrusters to actively position the liquid propellant away from 
tne lank vent. This method of pressure control was adequate for snort term missions 
anc oeemeo economically more feasible than tne weight penalty of additional insulation. 
(Ref. 7) Tne objective of tne present study is to predict the pressure response of a 
saturated llquio-vapor system wnsn unaergoJng a venting or depressurization process in 
zero gravity at low vent rates. 

Fig. 1 is a scnematic of a typical test container, with the liquio vapor interface 
assuming a hemispherical shape in zero-gravity. Fig. 2 is a schematic of the proposed 
venting model. The )-v interface is assumed to be planar, but with tne surface area of 
tne hemispherical interface, and the contents of Pie container are assumed to be at 
saturation conditions corresponding to Pv prior to venting, t<0. Upon initiation of 
venting, t>D, all properties are considered specially uniform but time dependent, except 
for tne liquid, whose temperature varies specially one-dimensionally as well. The 
interfacial temperature is tne saturation temperature corresponding to tne system 
pressure Pv. Tne analysis consists of applying tne appropriate governing equations to 
three control volumes; tne vapor, the liquid-vapor interface, and the liquid. Figures 3-5 



are schematics of these three control volumes. The vapor is treated as a lumped or 
uniform property control volume, and the conservation of mass and energy are applied, 
the interfacial mass transfer is found by applying the conservator of energy to the 
liquid-vapor interface. The liquid is treated as a semi-infinite planar solid in order to 
calculate the temperature gradient of the liquid at the interface. 

For purposes of comparison, an adiabatic model, which assumes no interfacial mass 
transfer, is constructed. The analysis, presented in Appendix O, is otherwise identical 
to that developed below. This model, when compared with the interfacial mass transfer 
model, will aid in evaluating the impact of interfacial mass transfer on the pressure 
response of the system. 

The pressure responses determined with the interfacial mass transfer and adiabatic 
models are compared with the results from previous models and with the experimental 
results obtained from the short duration drop tower tests conducted at the Lewis 
2 ero-gravity facility. 


NOMENCLATURE 


a thermal diffusivity, m2,sec 

a area, m2 

Cd discharge coefficient 

Cv specific heat at constant volume, J/kg-K 
S penetration depth, m 

h specific enthalpy, J/kg 

hfg neat of vaporization. J/kg 

k thermal conductivity, w/m-K 

m mass, kg 

n unit normal vector 

P pressure, N/m2 

q neat flux, w/m2 

R gas constant, m-N/kg-K 

T temperature, k 

t time, sec 

u internal energy, J 

u specific internal energy, J/kg 

v velocity, m/sec 

Subscripts: 

e vented vapor JJJHQINal Paqf to 

Of Poor, Quj*® 



1 liquid- vapor Interface 

1 liquid 

o Initial 

sat saturated conditions 

v ullage vapor 
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ANALYSIS 

The integral form of the continuity and energy equations for a control volume are 

used. 

The continuity equation is 

J dY 4 f fS-ndA = o (1) 

* A 

The volume v may oe assumed constant, since the actual volume changes due 
evaporation are small. Then, Eq. (l) Becomes 

^ ( 2 ) 

For tne vapor region, Eq. (2) becomes 

dmv/dt-mi-me (3) 

where mi is tne rate of generation of vapor at tne liquid-vapor interface, and me is the 
mass flow rate of tne vapor vented. For the liquid region 

omi/ot— ml (a) 


Tne energy equation is 


_d 

<k 



4 


j f kv • n M 

A 



O ( 5 ) 
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For purposes of the present analysis, it will oe assumed that: 

1. Heat transfer from the walls is negligibly small. 

2. No heat transfer takes place between the vapor and the i-v interface. 

3. The Internal energy in the vapor is specially uniform, varying only with time. 

4. The vapor volume 1$ constant (volume increases due to evaporation are 

neglected). 

5. The Interface surface area remains constant 

6. The liquid mass is large compared to the amount evaporated. 

7. All vapor properties are uniform at the state defined by Tv ana Pv. 

8. The interface temperature Ti-Tsat © Pv. 

9. Tne liquid-vapor mixture is initially saturated at Tv-Ti-Tsat ® Pv. 

For tne relatively short test times being modeled, along with the low venting rates 
assumed, these assumptions are reasonable. For longer test times, conduction from the 
walls must be taken into consideration. For the vapor t hen, Eq. (5) reduces to 

(Mv Uv ) 4 I^ c ky - - o ® 


Expanding Eq. (6): 

U v + /r»eky - h l » O 


now, assuming Cv-constant over a small temperature range, and substituting Eq. (3) into 
Eq. at 




Ar v 


+ ( Wv " K i. ) +* me Qta v - Mi/ ^ s O 


( 6 ) 


Expressions for ml ana me win now be developed. 

Tne mass flow rate tnrougn the vent, me, is determined by using a classical choked 
flow analysis (Ref. 9). Since the gas is venteo directly to a vacuum, the choked flow 
assumption is valid and the exiting mass flow rate is a function of upstream vapor 
properties only, given by: 


. Pv C. y Ar 


( 9 ) 


where Cd is an experimentally determined discharge coefficient and: 
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The rate of vapor generation, mi, is determined from the conservation of energy 
equation (Eq. (5)) applied to the liquid-vapor interface. Assuming no heat transfer to the 
vapor, all energy transferee to the interface by conduction in the liquid results In 
vaporizaion of liquid at the interface. Eq. (5) reduces to: 

* r*l iif j, (ID 

For relatively short periods, where the temperature boundary layer is small compared to 
any radii of curvature present at the Interface, the liquid may be treated as a semi-infinite 
planar solid. The surface area term, Al, will be the surface area of the hemisphere, the 
shape the interface takes in 2ero gravity. Referring to Fig. 5, the one dimensional 
conduction equation is: 

/«* A; ( i£)|^ o (12) 

Comolning equations (ll) and (12) gives 



Thus, the problem of determining to interfacial mass transfer is reduced to determine. g the 
temperature gradient of the liquid at the interface, which requires that the transient 
temperature distribution in the liquid near the 1-v interface be determined. If the liquid 
near the 1-v Interface can be considered to approximate a one-dimensional semi-infinite 
solid in lt*s thermal behavior the analytic solution for a step change in surface 
temperature, in connection with the finite form of Dunammel's superposition integral, can 
be used to determine the transient temperature distribution in the liquid. The time varying 
interface temperature is taken as the saturation temperature corresponding to the 
Instantaneous system pressure, which must be determined appropriately from the system of 
governing equations. 

Accordingly, the differential form of tne governing equation and the initial and 
oounaary conditions for the one-dimensional semi-infinite solid, initially at uniform 
temperature To and with a step change in surface temperature to Ti are: 


ST „ Z’T 

at T * 57* m 


r(*,o) - T* 
/ (0,1) * T l 


(15) 

(16) 
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rc-vo* t; 


(17) 


Tne solution is (Ref. 2> 


Tut') - Tj 

t. - n 


«r-f 


( 



(18) 


The Interface temperature, being tne saturation temperature corresponding to tne ullage 
pressure, will pe time varying in tne present case since tne pressure will cnange as tne 
tank is vented. Tnis time varying boundary condition Tl(t) is incorporated into tne solution 
using Dunammel*s superpostlon intregral (Ref. 2) in tne form: 


t 

o 

Here, 

- T 0 

- TiU) -T c 


(19) 


( 20 ) 


ano we let 


0Cx, ±) a 


O(Ki-) 

It) 


( 21 ) 


(i(x,t) is tne unsteady temperature resulting from a stepwise unit increase in surface 
temperature, relative to a uniform initial temperature, if tne increase is kept at zero until 
a certain time t-s, and at tnat instant raised to unity and maintained constant, tne new 
temperature ^(x,t) may be expressed in terms of p(x,t) as 



o t t< S 
j PU, ■*-*>, t>s 


( 22 ) 


me solution for ^(x,t) is given by Eq.(18), transformed to tne form of Eq. (21) as 
n*,<0 = S er-f c ( -Y~ pr ) (23) 

Solution of tne system of equations for tne venting problem will be performed in discrete 
time steps, and tne discrete form of Eq. (19) is given by: 

$• (X ,t) » (c) - + i: A ^ • (0 Of, i - so (24) 

/>>•/ * *' 



where 
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A * QiLSrr.') - P A ‘ ( $/»_, ) 


( 25 ) 


Here, n Is the total number of time steps into which the process has been divided, m is a 
running time index, i<m<n, and &&ijs the incremental change in surface temperature, 
related to the system vapor pressure. 

It is difficult to obtain a temperature gradient in the liquid at the interface to the 
desired degree of precision from the solution in the form of Eq. (24 )l Rather, the procedure 
followed here is to compute the instantaneous temperatures at a finite number of points in 
the liquid near the interface, using Eq. (24), and fit these points to a third order polynomial 
using a least squares fit. The polynomial used Is of the form: 


T - * A ♦ Bx f Cx 2 - + bx 3 ( 26 ) 


The temperature gradient of the liquid at the l-v interface, x-0, is there 


rr 

d. X 


X » o 


£ 


(27) 


The number and spacing of the nodes at which the temperatures of the liquid are to be 
calculated, and with which the coefficients A, B, C, and D in Eq. (26) will be determined, 
must next be specified Six nodes were taken arbitrarily as being sufficient to obtain the 
four coefficients in Eq. (26). Intuitively, nodes nearest to the l-v interface will give the 
most accurate value of the liquid temperature gradient at the l-v interface. The method 
used was to estimate a temperature penetration depth, , taken here to be the depth at 
which tne dimensionless temperature change computed by equation ( 18 ) is 95%. This is 
determined as: 


M5"» 

gr-f / 5 \ 

l 2 (a ±) iU ) 

(28) 

$- ' 

- 3 S -2(a 

(29) 


The actual penetration depth win be somewhat less than tnis value, since the actual system 
does not undergo a single step change in surface temperature, but rather a transient change 
in surface temperature. The six equally spaced nodes are taken to be within the 10% of this 
penetration depth nearest the l-v interface, shown schematically in Fig. 6. 

now tnat the temperature of the liquid at each of the six nodes near the l-v interface is 
known, the constants A ££, and D of Eq. (26) may be determined. A least squares 
algoritnm was used (Ref. 4), which determines tn* polynomial coefficients which minimize 
the error between the data points and the polynori.ial. 
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Appendix A describes a test program devised to evaluate the effect of the fraction of 
penetration depth used when fitting a polynomial by computing the accuracy of the 
polynomial in predicting the temperature gradient at the 1-v interface. The temperature 
gradient obtained with the above procedure is compared with the analytical value for a 
single step change in surface temperature,, being the most severe test possible. This is 
done for different fractions of the penetration depth. Figures Ai and A2 show that with 
the nodes spaced in a region of 10% of the penetration depth from the surface and using a 
third order order polynomial, an error of less than 0.5% in temperature gradient at the 
surface is obtained. 

For the adiabatic model, the mass transfer at the interface is taken as zero, and the above 
analysis for the interfacial mass transfer is not used. 

When combined with the proper initial conditions, equations (3), (8), (9), and (13), along with 
the liquid temperature distribution, provide a complete description of the vapor space. 

These equations were numerically solved by computer. A program listing and description is 
included in appendix 6 . A comparison of these results with the experimental data 
available to date is presented below. 


RESULTS 

The model described above differs primarily in two respects from previous models used to 
predict the pressure response of an initially saturated liquid vapor mixture vented to a 
vacuum in zero gravity. The most significant difference is the procedure used to 
approximate tne interfacial mass transfer. The present model assumes the liquid to be a 
semi-infinite solid with a planar surface and a transient surface temperature determined 
from tne coupling between the liquid conduction process and the vapor behavior. 

Ounammel's superpostion integral is used to incorporate tne effect of a transient surface 
temperature in computing the liquid temperature profile. The interfacial mass flux is then 
determined from the temperature gradient at the liquid-vapor interface. 

The second difference from past models 1$ that the vapor temperature is not assumed to be 
at tne saturation temperature corresponding to the vapor pressure. This now couples the 
energy and continuity equations for the vapor system and makes for a more difficult 
numerical solution. Tne effect of this change in assumption can be seen in figures 7 and 8, 
wnere both the mean vapor temperature and tne instantaneous saturation temperatures are 
plotted for two test runs. The difference between tne vapor temperature and the saturation 
temperature can be as much as 30 degrees K. The vapor temperature Is higher than the 
saturation temperature and is thus superheated. Since me is Inversely proportional to vapor 
temperature, hlgner vapor temperatures result in slightly lower vent rates, and thus slower 
ullage pressure crop. 

Comparison between the pressure response predicted by the present model, the present 
adiabatic model, and previous models (Ref. 1) are given in Table l, together with 
measurements obtained previously (Ref. l). The data in Table l snows that the proposed 
model gives pressure responses closer to the experimental data than does any previous 
model. The data in Table l also shows that botn the present model and previous models 
incorporating interfacial mass transfer yield better results tnan does the adiabatic model, 
which assumes no interfacial mass transfer. It is evident that interfacial mass transfer 
must be considered when using low vent rates such as the ones used in this study. Hence, 
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it may be concluded that the proposed model better approximates interfacial mass transfer 
than previous models, but the sizeable error when compared to the experimental data 
indicates that certain elements are still lacking in the description of the process. It Is also 
possible that the experiments themselves should reexamined. 

Additional detailed transient behavior of Runs 2 and 4 in Table l are plotted in figures 
9-14, with system pressures in Figs. 9 and 10, vent rates in Figs. 11 and 12, and 
evaporation rates in Figs. 13 and 14. Run 4 has a discharge area 2.2 times that for Run 2, 
approximately the same initial volume, and an initial pressure approximately 10% higher. 

This is consistent with the higher pressure drop rate, higher vent rate, and higher 
evaporation rate that occurs wlith Run 4. 

Evaluation of this model assumes that the experimental data accurately describes the 
system being modeled. The small test vessels used would tend to make the geometry of the 
system Important. The flow coefficients, Cd, were experimentally determined, and there is 
no way of evaluating their accuracy. Future experiments should be conducted before 
making a final evaluation of the model proposed here. 


CONCLUSION 


An analytical model was constructed to predict the pressure response of cylindrical 
containers initially filled witn a saturated liquid-vapor mixture vented to a vacuum under 
2ero gravity conditions. The response predicted by this model was compared to that of 
previous models and to the experimental data obtained at the NASA Lewis Researach 
Center. 

Previous models predicted too large a pressure drop. The model proposed here gives a 
pressure response closer to the experimental data than other models, but still predicts too 
large a pressure d'op. This means that the present model still underestimates the amount 
of interfacial mass transfer. Higher rates of evaporation will yield a lower pressure drop 
in the system. An additonal source of vapor formation not considered in the present model 
is the thin liquid layer existing at the liquid-vapor-solid triple interline formed by a 
hemispherical liquid-vapor interface. It can be expected that rapid evaporation would take 
place in tnis region, involving conduction effects from the container walls (neglected in the 
present analysis). This would reduce the pressure drop predicted by the model, with 
pemaps better agreement with experiments conducted to date. 

Future experiments might be considered for comparison with the present model in 
which the presence of the triple interline would be minimized by using larger size vessels 
and by . conducting the experiments at standard earth gravity. 



TABLE 1. SUnrtARY OF PARAHETERS AND RESULTS 



Test Initial 
Run vapor 
No voluno 

Norzla 
(liana tar 

Discharge 

coeffic. 

Initial 

ullage 

pressure 

Initial 

ullage 

tenp. 

Final 

exo. 

pressure 

Final 

analy. 

press. 

Final 

past 

analy. 

press. 

Final Dinon- Dinen- 

adiabatic sionless sionless 
press. anal. exp. 

press drop press dr-' 

nJ 

n 

Cd 

kPa 

K 

kPa 

kPa 

kPa 

kPa 



1 1.93E-4 

0.406E-3 

0.64 

89.6 

294.3 

86.2 

85.2 

81.6 

83.2 

0.07 

0.06 

2 2.01 

0.889 

0.69 

87.9 

294.7 

70.3 

64.4 

56.3 

56.1 

0.31 

0.25 

3 1.90 

1.07 

0.86 

91.0 

293.7 

60.7 

46.8 

40.7 

33.6 

0.48 

0.33 

4 1.93 

1.32 

0.875 

97.2 

296.5 

53.8 

37.9 

29.4 

21.8 

0.62 

0.46 

5 1.93 

1.93 

0.77 

101.0 

295.4 

41.4 

21.4 

13.1 

5.3 

0.78 

0.57 
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LIQUID-VAPOR INTERFACE 



VAPOR 

LIQUID 

NODES 



Figure 6. - Location of nodes in liquid. 



Pun No. 2 
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Temp., K T vapor 






Pressure, kPa 


Run 2. 
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Time, secs. 

Figure 9. - Transient System Pressure Response 







Run No. 2 
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Evaporation rate, g/sec 

A M 



-I 


Time, secs. 

Figure 13. - Transient interface Evaporation Rate 


Run No. 4 


Evaporation rate, g/sec 



Time, secs 

Figure 14. - Transient Interface Evaporation Rate 
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EXAMINATION OF THIRD ORDER CURVE FIT ACCURACY 


A test was run to determine the accuracy of the third order least squares curve 
fit used in calculating the interfacial mass transfer. The test also determined the 
spacing of the nodes in the liquid which would give the best curve fit. The 
temperature distritxtion in a semi-infinite solid with constant surface temperature is 
derived in the ANALYSIS and given by Eq. (18): 


TiKi'i - Tc 

To - n 



m 


From this the temperature gradient at x=0 is- 

ctT“ I _ ! c - / c (A2) 

tix /x=0 (wort y ; *. 


Equation (A2) represents an exact solution for the gradient. An approximate solution 
is Detained via the method described in ANALYSIS. The penetration depth is 
calculated. A percentage of this depth near the surface is then divided into six 
equally spaced nodes at which the temperature is calculated. First, second, and third 
order curves are fit to the data obtained using different percentages of the 
penetration depth. As can be expected, the calculated gradient and intercepts were 
most accurate when the nodes were space closest to the interface, i.e. a small 
percentage of the penetration depth. Figures A1 and A2 show that using a third 
order polynomial with nodes very close to the interface give the best gradient and 
intercept results. 

Note that equation (A2) gives the exact gradient for a semi-infinite solid with a step 
change in surface temperature. This can not be used in determining interfacial mess 
transfer in the proposed model of this report since the surface temperature in reality 
is a f unction of time. 







ORIGINAL PAGE IS 
APPENDIX B OF POOR QUALITY 

COMPUTER ALGORITHM 


The conputer algorithm used to numerically solve the governing equations 
consists of a main program and eight subroutines. The basic outline of the 
numerical solution is as follows. At time= t, the vapor temperature and mass( 
y(i), y( 2) ) are known, along with the ullage volume, which is constant. Thus 
the state of the vapor aid the interface are defined and all thermodynamic 
properties can be determined. With the state at time=t conpletely defined, 
values of y(i) and Y(2) at t*t ♦ 0.05 are found by solving the governing 
differential equations by a fourth-order Runge-Kutta method, with the values of 
y(i) and Y(2) now determined at time«t ♦ 0.05, this state is now completely 
defined, and the algorithm can be incremented by one time step and repeated. The 
following is a brief description of the function of each subroutine. 

RUNSE - A fourth order Runge-Kutta algorithm to solve first order differential 
equations with non-constant, coefficients. This routine uses a fixed time step, 
with the time step being the independent variable. 

DERY - calculates the derivatives of y(i) and Y(2) with respect to time for use 
in the RUNGE algorithm. 

PROPS - Determines the necessaary thermodynamic properties of tne working fluid, 
given vapor temperature, mass, and volume. The four basic equations used to 
calculate the properties are; vapor-pressure equation, equation of state, 
density of saturated liquid, and neat capacity of vapor(Ref.3 ). All properties 
can be determined from these equatlons(App.C). 

NEWTTS - The vapor-pressure equation is of the form P=f(Tsat). This routine 
uses the Newton-Rap son metnod(Ref .4) to solve this equation for Tsat, given P. 

newtv - The equation of state is of the form P«f(v,Tv). This routine uses the 
Newton-Rapson method to solve the equation of state for the specific volume v, 
given P and Tv. These values of v are needed in PROPS to calculate internal 
energy and enthalpy. 

mass - Determines the rate of mass transfer across the liquid-vapor interface. 

As discussed in analysis, the liquid temperature gradient at the interface Is 
needed to compute the interfacial mass transfer. Duhammel's superposition 
integral and the one dimensional conduction equation for a semi -infinite solid 
with a step change in surface temperature are used to compute the temperature of 
i>'- liquid at various depths near the interface. A third order least squares 
curve fit(Ref.A) is used to find the best curve through these points and thus 
the surface temperature gradient. 

slud - Along with SUR, solves the system of equations describing the third 
order least squares curve fit. This routine computes the LU decomposition of the 
coefficient matrix. 

SLIR - Computes the solution to the system of linear equations AX*B using 
iterative refinement. SLUD and SLIR are called from tne MTS Numerical Analysis 
Library(Ref .5). Similar routines are readily available for users not on the MTS 
network (Ref .4). 



FORTRAN 

Symbol 

Engineering 

Symbol 

Description 

Units 

AS 

As 

Interface surface area 

ft2 

AT 

At 

Nozzle cross sectional area 

ft2 

CO 

Cd 

Discharge coefficient 

- 

CVTVP 

cv 

Specific heat of vapor § TV,P 

ft-lbf/slug-R 

HFGTS 

hfg 

Enthalpy of evaporation 8 TS 

ft-lbf/slug 

HVTSP 

hg 

Enthalpy of vapor 9 TS,P 

ft-lbf/slug 

HvTVP 

h 

Enthalpy of vapor 8 TV,P 

ft-lbf/slug 

KLTS 

k 

Thermal conductivity of liquid 8 TS 

lbf/sec-R 

re 

me 

Macs flow rate of vapor vented 

slug/sec 

MI 

mi 

Mass flux across 1-v interface 

slug/sec 

P 

P 

Ullage pressure 

lbf/in2 

PR 

Pref 

Peference pressure 

lbf/in2 

R 

R 

Ideal gas constant 

psi-ft3/lbm-R 

T 

t 

Time 

seconds 

TC 

Tc 

Critical temperature 

R 

TR 

Tref 

Reference pressure 

R 

TS 

Tsat 

Saturation temperature 8 P 

R 

UVTRPR 

uref 

Reference internal energy f TR, PR 

ft-lbf/slug 

U'/TSP 

ug 

Internal energy of vapor 8 TS,P 

ft-lbf/slug 

UVTVP 

u 

Internal energy of vapor 8 TV,P 

ft-lbf/slug 

VU 

Vu 

Ullage volume 

ft3 

VLTSP 

vf 

Specific volume of liquid 8 TS,P 

ft3/lbm 

WTSP 

vg 

Specific volume of vapor 8 TS,P 

ft3/lbm 

WTVP 

V 

Specific volume of vapor 8 TV,P 

ft3/lbm 

WTVPR 

V 

Specific volume of vapor 8 TV, PR 

ft3/lbm 

Y(l) 

Tv/ 

Temperature of ullage vapor 

R 

Y(2) 

mv 

Mass of ullage vapor 

lbm 

YP(1) 

dTv/dt 

Time rate of change of vapor temp. 

R/sec 

YP(2) 

cfrnv/dt 

Time rate of change of vapor mass 

lbm/sec 
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Listing of MAIN*... at 10:07:26 on APR 6, 1984 for CCid«SS3X 


1 

2 

3 

4 

5 

6 

7 

8 
9 


C 

C 

C 

C 

C 




READ IN CONSTANTS FOR FREON- 11 VAPOR PRESSURE CURVE, 
OF STATE, HEAT CAPACITY OF THE VAPOR, AND DENSITY OF 
LIQUID. 

BLOCK DATA 

REAL A(6),B(6),C(6),D(6),E(6),F(6),R,CK, 

1 SB , P , MS , TCRI T , TR , PR , CD , CC , AT , UV , AS , KLTS 


EQUATION 

SATURATED 


10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 


COMMON/ALPHA/A, B,C,D,E,F,R,CC, SB, CD, AT, UV,TCRIT,TR, PR, 

1 UVTRPR , AS 

C 

DATA A/0.0,-3. 126759,-0.025341 ,0.00 1687277, -2. 35893E-5, 

1 1 . 057504E8/ 

DATA B/0. 0,0. 00 13 18523, 4.875 12 IE-5,- 1.805062E-6,2.448303E-8, 
1 -9.472103E4/ 

C 

DATA C/0.0,-35.76999, 1.220367, 0. 0,-1. 478379E-4, 0.0/ 

C 

DATA D/42.14702865,-4344.343807,-12.84596753, 

1 0.004008372507,0.0313605356,862.07/ 

C 

DATA E/34.57,57.63811,43.63220,-42.82356,36.70663,0.0/ 

C 

DATA F/0 .023815,-336.80703,2 . 798823E-4 , -2 . 123734E-7 , 

1 5. 9990 18E- 11,0.0/ 

C 

DATA R,TCRIT,SB,CC/0. 0781 17,848.07,0.0019,-4.5/ 

DATA TR , PR , UVTRPR/4 19.67,0.74137,2032163.0/ 

END 

C 

C BEGIN MAIN PROGRAM. 

C 

C 

REAL Y(2),YP(2),A(6),B(6),C(6),D(6),E(6),F(6),R,CK, 

1 SB, P, MS, ERF ( 150,2) ,TCRIT,TR,PR,CD,CC,AT,UV,AS,KLTS 
REAL 2(4) 

REAL MI ,ME,WME,WMI ,WTS,TS,WMASS 
C 

COMMON/ALPHA/A, B,C,D,E,F,R,CC, SB, CD, AT, UV,TCRIT,TR, PR, 

1 UVTRPR, AS 
C 

C READ IN ERROR FUNCTION VALUES FOR USE IN SUBROUTINE MASS. DATA 
C LOCATED IN FILE 'ERF'. 

C 

DO 22 K= 1 , 102 

READ (7, 34) ERF ( K , 1 ) , ERF (K , 2 ) 

34 FORMAT ( 2F20. 9) 

22 CONTINUE 

C 
C 
C 

C TR, PR, UVTRPR INITIALIZED IN SUBROUTINE PROPS 
C 

C INITIALIZE T,Y( 1 ) ,Y{2) , AT T=0.0 SECONDS. 

C 

C SET AT AND CD, THE VARIABLES WHICH CONTROL VENT FLOW RATE. 

C SET UV, THE ULLAGE VOLUME 



Listing of MAIN+... 


at 10:07:26 on APR 6, 1964 for CCid«SS3X 


59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
61 
82 

83 

84 

85 

86 
87 
86 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
1 1 1 
112 

113 

114 

115 

116 


C 

C UNITS: T IN SECS.,Y(1) IN RANRINE, Y(2) IN LBM, AT IN FT2, 

C UV IN PT3. 

C INITIALIZE VARIABLES AT T-0.0 
C 

T<=0.0 

Y( 1 )=531 . 74 

Y( 2 )*»0 . 00237026 ORIGINAL PAGE (S 

AT-0. 0000314902 OF POOR QUALITY 

UV=0. 0068 15371 
CD=0. 77 
C 

C WRITE OUT INPUT VALUES 

WRITE (6,38) Y ( 1 ) ,'Y ( 2 ) , AT , UV , CD 

38 FORMAT (’ INITIAL VAPOR TEMPERATURE IS’ ,2X,F8.3, 'RANKINE'/ 

1 'INITIAL VAPOR MASS I S ’ , 2X,E1 1 . 5 , 2X, ' LBM’ / 

1 'NOZZLE AREA IS ' , 2X,E1 1 . 5 , 2X, ' SQUARE FEET’/ 

1 'ULLAGE VOLUME IS ’ , 2X, E 11 . 5 , 2X , ' CUBIC FEET'/ 

1 'DISCHARGE COEFFICIENT IS',2X,F6.2) 

C 

C WRITE OUT HEADINGS 
C 

WRITE(6, 71 ) 

WRITE(6 , 72 ) 

WRI TE ( 6 , 7 3 ) 

73 FORMAT ( ' ') 

71 FORMAT( 'TIME' ,2X, 'T VAPOR’ ,4X, 'TSAT' ,5X, 'P VAPOR' ,3X, 

1 'VAPOR MASS’ ,2X, 'VENT RATE' ,3X,.'EVAP RATE' ) 

72 FORMAT( 'SECS' ,3X, 'KELVIN' , 4X, ’ KELVIN ’ , 3X, ’ PASCALS ’ ,6X,'KG' , 

1 8X, ' KG/SEC ' ,BX, 'KG/SEC' ) 

! 

c 

C USING A FOURTH ORDER RUNGE KUTTA METHOD TO EVALUATE THE INTEGRALS 
C THE FOLLOWING LOOP WILL BE RUN THROUGH 60 TIMES WITH A TIME STEP : ■ 
C OF 0.05 SECONDS. TOTAL TEST TIME BEING 3.0 SECONDS. ;J 

C if 

DO 23 KL® 1 , 60 

CALL RUNGE ( Y , T , YP , P , MI , ME , ERF , TS ) 

C 

C CONVERT UNITS FROM ENGLISH TO MKS AND WRITE OUT RESULTS 
C 

WTEMP® (Y(1)-459.67) *5/9“ 17.77778+ 273.14 
WTS® (TS-459. 67 ) *5/9 -17.77778 ♦ 273.14 
WMASS«Y(2)/2.205 
WP=P*6895.0 
WMI®MI* 14 .59 
WME=ME* 14.59 
C 

C INCREMENT TIME 
C 

T=T+0 . 05 

WRITE( 6,24) T,WTEMP,WTS,WP,WMASS,WME,WMI 
24 FORMAT (F4 ,2,2X,F8.4,2X,F8.4,2X,F8.2,2X,F10.8, 1X,E 1 0 . 4 , IX, El 1 . 

23 CONTINUE 

STOP 
END 
C 

C SUBROUTINE RUNGE RUNS A 4TH ORDER RUNGE-KUTTA METHOD TO NUMERI- , 
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at 10:07:26 on APR 6, 1984 for CCid*SS3X 


117 

118 
H9 
120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 
167 
166 

169 

170 

171 

172 

173 

174 


C CALLS SOLVE THE DIFFERENTIAL EQUATIONS GOVERNING THE SYSTEM. 
C 

SUBROUTINE RUNGE(Y, T,YP, P,MI ,ME,ERF,TS) 

REAL A{6),B(6),C(6),D(6),E(6),F(6),R,CK, 

1 SB , P , MS , ERF (150,2) , TCRI T , TR , PR , CD , CC , AT , UV , AS , KLTS 
REAL Z (4) 

C 

COMMON/ALPHA/A, B, C, D,£,F,R,CC, SB, CD, AT, UV, TCRIT.TR, PR, 

1 UVTRPR , AS 

REAL ENDRKS, MI ,ME 

REAL K0(2),K1(2),K2(2),K3(2),Y(2),YP(2) ,NEWY(2) 

H=0. 05 

C TIME STEP, H, SET AT 0.05 SEC 
C 

C COMPUTE FIRST APPROX OF SLOPE 
C 

ENDRKS = 1 . 0 

CALL DERY ( Y , T , YP , P , MI , ME , ENDRKS , ERF , TS ) 

DO 90 J=1 ,2 

K0 ( J )=H*YP( J ) 

90 CONTINUE 

C 

C SECOND APPROX OF SLOPE 
C 

ENDRKS =0.0 
Z(2)=Y(2)+K0(2)/2. 

Z( 1 )=Y( 1 )+K0{ 1 )/2.0 
V=T+H/2 . 0 

CALL DERY { Z , V , YP , P , MI , ME , ENDRKS , ERF , TS ) 

K 1 ( 1 )=H*YP{ 1 ) 

K1.(2)=H*YP(2) 

C 

C THIRD APPROX OF SLOPE 
C 

Z ( 1 ) =Y( 1 ) +K1 ( 1 )/2 . 0 
Z ( 2 ) = Y ( 2 ) +K 1 ( 2 ) /2 . 0 

CALL DERY ( Z , V , YP , P , MI , ME , ENDRKS , ERF , TS ) 

K2( 1 )=YP( 1 )*H 
K2(2)=YP(2)*H 
C 

C FOURTH APPROX OF SLOPE 
C 

Z( 1 )*Y{ 1 )+K2 ( 1 ) 

Z(2)=Y(2)+K2(2) 

V=T+H 

CALL DERYU,V,YP,P,MI , ME, ENDRKS, F,TS) 

K3( 1 )=H*YP( 1 ) 

K3(2)=H*YP(2) 

C 

C PREDICT FUTURE Y BASED ON AN AVERAGE SLOPE 
C 

DO 93 M«1 ,2 

Y(M)=Y(M)+(K0(M)+2*K1(M)+2*K2(M)+K3(M) )/6.0 
NEWY (M)*Y (M) 
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C 

C 


93 


CONTINUE 



Listing of MAIN'*-. . . 


at 10:07:26 on APR 6, 1984 for CCid*SS3X 


175 

17fr 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 
153 

194 

195 

196 

197 
196 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 


RETURN 

end ORIGINAL PAGE 19 

c OF POOR QUALITY 

C SUBROUTINE EVALUATES THE VALUES dY(1)/dt AND dY(2)/dt FOR EACH 
C CALL FROM THE SUBROUTINE RUNGS. 

C 

SUBROUTINE DERY( Y,T, YP,P,MI , ME , ENDRKS , ERF , TS ) 

REAL ENDRKS 

REAL A(6) ,B(6) ,C(6) r D(6) ,E(6) , F(6) ,R,CK, 

1 SB, P, MS, ERF ( 150,2) , TC RI T , TR , PR , CD , CC , AT , UV , AS , KLTS 
C 

COMMON/ALPHA/A, B,C,D,E,F,R,CC, SB, CD, AT, UV,TCRIT,TR, PR, 

1 UVTRPR , AS 

C 

REAL Y ( 3 ) , T , YP ( 3 ) , MI , ME 
C 

C CALL SUBROUTINE PROPS TO FIND THERMODYNAMIC PROPERTIES OF THE 
C LIQUID AND VAPOR, GIVEN Y(1),AND Y(2); THE TEMPERATURE AND MASS C 
C THE VAPOR. 

C 

CALL PROPS ( Y , T , YP , P , ME , HVTVP , UVTVP , HVTSP , HFGTS , TS , CVTVP , KL1 
C 

C CALL SUBROUTINE MASS TO COMPUTE THE MASS FLOW RATE ACROSS THE 
C LIQUID-VAPOR INTERFACE 
C 

CALL MASS ( Y , T , TS , HFGTS , KLTS , VLTSP , MI , ENDRKS , ERF ) 

C 

C COMPUTE DY(1)/DT AND DY(2)/DT, THE DERIVATIVES OF VAPOR TEMPERATl 
C AND VAPOR MASS WITH RESPECT TO TIME. j 

C 

YP ( 1 ) = ( HVTSP-UVTVP ) *MI / ( Y ( 2 ) ♦CVTVP ) + ( UVTVP-HVTVP ) *ME/ ( Y ( 2 ) *C\| 
YP(2) = (MI -M£)*3'2. 174 
C 
C 


RETURN 

END 


212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 


C 

C SUBROUTINE PROPS COMPUTES THE THERMODYNAMIC PROPERTIES OF THE 
C WORKING FLUID, GIVEN THE VAPOR TEMPERATURE AND MAS". 

SUBROUTI NE PROPS ( Y , T , YP , P , ME , HVTVP , UVTVP , HVT r ,P , HFGTS , 

1 TS, CVTVP, KLTS, VLTSP) 

C 

COMMON/ALPHA/A, B,C,D,E,F,R,CC, SB, CD, AT, UV,TCRIT,TR, PR, 

1 UVTRPR AS 

REAL Y(2)'yP(2),A(6),B(6),C(6),D(6),E(6),F(6) ,ERF( 150,2) 

REAL XV( 4 ) ,XT( 4 ) , WV( 4 ) ,WT(4 ) 

REAL KLTS , ME , XCV ( 4 ' 

C 

C CRITICAL TEMPERATURE AND RELATIVE TEMPERATURE AND PESSURE OF FRECj 

C 1 

*\ 8.07 

P .74317 
7k -427,0 
C 

C COMPUTE SPECIFIC VOLUME OF ULLAGE VAPOR, FT3/LBM. 

C 


WTVP»UV/Y(2) 



at 10:07:26 on APR 6, 1984 for CCid=SS3X 


Listing of MAIN*. . . 


233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 

274 

275 

276 

277 

278 

279 

21m 

281 

282 

283 

284 

285 

286 
287 
286 

289 

290 


C 

C COMPOTE ULLAGE PRESSURE FROM EQUATION OF STATE, KNOWING TEMPERA Tl 
C OF ULLAGE AND SPECIFIC VOLUME OF ULLAGE VAPOR , UNITS OF P ARE PS 
C 

P=R*Y(1)/(WTVP-SB) ■*-(A(2)+B(2)*Y( 1 )+C(2)*EXP(CC*Y ( D/TCRIT 
1 /( (WTVP-SB)**2) +(A(3)+B(3)*Y(1)+C(3)*EXP(CC*Y(1)/TCRIT 

1 /( (WTVP-SB)**3) +(A(4)+B(4)*Y(1))/((WTVP“SB)**4) 

1 ♦ (A(5)+B(5)*Y(1)+C(5)*EXP(CC*Y(1)/TCRIT))/((WTVP-SB)**5) 

C 

3 FORMAT (3F It. 9) 

C CALL THE NEWTONS METHOD SUBROUTINES TO FIND TS, WTSP, WTVPR, AND 
C WTSPR. THESE SPECIFIC VOLUMES ARE NEEDED FOR THE CALCULATION 
C OF ENTHALPY AND INTERNAL ENERGY. 

C NEWTTS USES NEWTONS METHOD TO SOLVE THE VAPOR PRESSURE EQUATION 
C FOR TSAT, GIVEN PSAT. NEWTV SOLVES THE EQUATION OF STATE FOR 
C SPECIFIC VOLUME, GIVEN TEMPERATURE AND PRESSURE OF THE VAPOR. 

C UNITS ARE: TS IN DEGREES RANRINE, SPEC. VOL. IN FT3/LBM 
C 

VTOL=0. 005 
TSTOL=0.5 

CALL NEWTTS (TSTOL , P , TS ) 

CALL NEWTV (VTOL,P,TS, WTSP) 

CALL NEWTV (VTOL, PR, Y ( 1 ) , WTVPR) 

CALL NEWTV ( VTOL, PR ,TS, WTSPR) 

C 

C ASSIGN TEMPORARY VALUES TO SPECIFIC VOLUMES AND TEMPERATURES TO CC 
C INTERNAL ENERGY AND EHTHALPY 
C 

XT(1)=TR 
XT ( 2 ) =TS 
ZT( 3 )=Y ( 1 ) 

XT ( 4 ) =TS 
XV ( 1 )= WTSPR 
XV ( 2 ) =WTVP 
XV( 3)=WTSP 
XV ( 4 ) = WTVPR 
C 

C EVALUATE THE INTERNAL ENERGY INTEGRALS IN THIS LOOP 
C UNITS OF WV AND WX ARE FT-LBF/SLUG OR FT2/SEC2 
C 

DO 6 1=1,4 
I F ( I . LE . 2 ) TEMP=TS 
I F ( I . GT . 2 ) TEMP=Y{ 1) 

BETA=TEMP*CC/TCRIT 

WV(I)=((A(2)-C(2)*( BETA- 1 . ) * EXP ( BETA) ) /(XV( I ) -SB) 

1 +(A(3)~C(3)*( BETA- 1 . ) *EXP ( BETA) )/( 2. * (XV( I )-SB) **2 . ) 

1 +A(4)/(3.*(XV(I)-SB)**3.) + (A(5)~C(5)*( BETA” 1 . ) *EXP ( BET* 

1 )/(4.*(XV(I )-SB)**4. ) )* 144.0*32. 174 

C 

C 

WT(I ) = (F(1)*XT(I )“F(2)/XT(T N 7 1 *:.T( I ) **2 . )/2 . 

1 ♦(F(4)*XT(I)**3. )/3. ♦(F4-- , -,.)/4.)*778. 17*32.174 

C 

C 

6 CONTINUE ‘ 


ORIGINAL PAGE 12 
OF POOR QUALITY 



Listing of MAIN+... at 10:07:26 on APR 


6 , 


1984 for CCid*SS3X 


ORIGIN At. PAGE 18 i 
OF POOR QUALITY 


291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 


C 

C COMPUTE INTERNAL ENERGY AND ENTHALPY USING THE VALUES OP WV AND 1 
C UNITS ARE PT2/SEC2 OR FT-LBF/SLUG 
C 

UVTVP-WV(2)-WV(4)+WT(3)'WT( 1 )+UVTRPR 
UVTSP=WV ( 3 ) -WV ( 1 ) +WT ( 2 ) -WT ( 1 ) +UVTRPR 
HVTVP=UVTVP + P* 144. 0*32. 174*WTVP 
HVTSP*UVTSP + P*144.0*32. 174*WTSP 
C 

C COMPUTE SPECIFIC VOLUME OF LIQUID IN FT 3 /SLUG 
C 

CON* 1- ( TS/TCRIT) 

RHOL*E( 1 ) +E(2)*CON**(1./3.) +E(3)*CON**(2./3. ) + 

1 E(4)*CON +E(5)*CON**(4./3. ) 

VLTSP=32. 174/RHOL 
C 

C COMPUTE DP/DT 

C 

C 

DPDT* ( -D ( 2 ) * ALOC ( 10.0)/(TS**2. ) ♦D(4)*ALOG( 10. ) +D(3)/TS 
1 -D(5)*D(6)*ALOG(D(6)-TS)/(TS**2. ) +D( 5)/TS) *EXP(ALOG( 10. ) 

1 *(D(1)+D(2)/TS +D(4)*TS) +D(3)*ALOG(TS) +D(5)* (D(6)-TS)* 

1 ALOG(D( 6 )-TS)/TS) 

C 

C COMPUTE ENTHALPY OF FORMATION 
C 

C IN FT2/SEC2 
C 

HFGTS=TS* ( WTSP- (VLTSP/32 . 1 74 ) ) *DPDT* 144 . 0*32 . 1 74 
C 

C COMPUTE K, THERMAL CONDUCTIVITY OF THE LIQUID 
C UNITS ARE LBF/SEC-R. A LINEAR CURVE FIT IS USED. 

C 

KLTS= ( 0 . 1 1 1 562“TS*0 . 000 1 15)*0. 216158 
C 

C COMPUTE CV, THE SPECIFIC HEAT, AT TEMPERATURE OF THE VAPOR 
C 

CV0*F( 1)+F(2)/(Y< 1)**2.) +F(3 ) *Y( 1 ) +F(4 )*Y( 1 )**2.+F(5)*(Y( 
C 

C THE DO LOOP EVALUATES AN INTEGRAL TO FIND SPECIFIC HEAT AT TV 
C RELATIVE TO THE SPECIFIC HEAT AT T~ RELATIVE 
C 

DO 357 L=2,4,2 

XCV(L )*Y ( 1 ) * ( “CC/TC ) **2 . *EXP(CC*Y ( 1 ) /TC ) * ( -C( 2 ) /(XV(L)“SB) 
1 -C(3)/(2.*(XV(L)-SB)**2. ) -C ( 5 ) /(4* (XV(L)“SB) **4 . ) )* 144 . *3 

357 CONTINUE 
C 

C XCV IN UNITS OF FT2/SEC2-R, CONVERT CV0,R TO THOSE UNITS 
C 

CV0=CV0*778. 16*32. 174 
ReR* 1 44 . *32. 174 
C 

C COMPUTE CV, THE SPECIFIC HEAT CONSTANT 
C 

CV*CV 0 + XCV ( 2 ) “ XCV ( 4 ) 

CVTVP=CV 

C 

C COMPUTE THE MASS FLOW RATE THROUGH THE N022LE BASED ON THE BULK 



Listing of MAIN*. . . 
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349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 


C 

C 

C 


C 

C 


PROPERTIES OF THE VAPOR 
IN UNITS OF SLUGS/SEC 
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1 


CP-R + CV 

ALPHA® (CP/CV)* *0.5 *(2/((CP/CV)+1.))**(((CP/CV)+1.)/(2. 
*((CP/CV)-1.))) 

ME=CD*ALPHA*P*AT* 144 ./( (R*Y( 1 ) )**0.5) 


R«R/( 144. *32. 174) 
RETURN 
END 


C 

C 

C 

SUBROUTINE NEWTV( ERROR, PRESS, TEMP, Z) 

REAL A(6),B(6) ,C(6) ,D(6) ,E(6) ,F(6),R,CR, 

1 SB, P, MS, ERF ( 150,2) ,TCRIT,TR,PR,CD,CC,AT,VU,AS,KLTS 
C 
C 

COMMON/ALPHA/A, B, C,D,E,F,R,CC, SB, CD, AT, VU,TC,TR, PR, 

1 UVTRPR , AS 

C 
C 
C 
C 

C THIS ROUTINE USES NEWTONS METHOD TO FIND THE ROOTS OF THE 
C EQUATION OF STATE EQUATION, THE SPECIFIC VOLUME. 

C 

C INITIAL GUESS FOR SPECIFIC VOLUME 
C 

X«R*TEMP/PRESS 

C 

C PERFORM NEWTONS METHOD UNTIL ERROR IS LESS THAN VTOL 
C 

DO 40 J= 1 , 7 

Z*R*TEMP/(X-SB) + ( A ( 2 ) +B ( 2 ) *TEMP+C ( 2 ) *EXP(CC*TEMP/TC) ) 

1 / ( (X-SB) **2 ) +(A(3)+B(3)*TEMP*C(3)*EXP(CC*TEMP/TC)) 

1 / ( (X-SB) **3 ) +(A(4)+B(4)*TEMP)/((X-SB)**4) * 

1 ( A( 5 ) + B( 5) *TEMP+C ( 5 ) *EXP(TEMP*CC/TC) )/( (X-SB) ** 5) “PRESS 
CON=CC*TEMP/TC 

DZDV=-(R*TEMP)/(X-SB)**2. -2 . * ( A( 2 )+B ( 2 ) *TEMP+C (2 ) *EXP(CON) ) 
1 /(X-SB)**3. -3.*(A(3)+B(3)*TEMP+C(3)*EXP(CON))/(X-SB)»*4. 

1 -4.0*{A(4)+B(4)*TEMP)/(X-SB)**5. -5 . * ( A( 5)+B( 5) *TEMP 

1 +C(5)*EXP(CON) )/(X-SB)**6. 

C 

C COMPUTE NEW SPECIFIC VOLUME 
C 

X=X-Z/DZDV 

IF (Z/DZDV.LT. ERROR) GO TO 40 
40 CONTINUE 

RETURN 
END 

SUBROUTINE NEWTTS (ERROR , PRESS ,X ) 

REAL A(6) ,B(6) ,C(6) ,D(6) ,E(6) ,F(6) ,R,CK, 

1 SB , P , MS ,ERF (150,2) ,TCRIT , TR , PR, CD,CC, AT, VU, AS ,KLTS 
C 
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407 

408 

409 

410 

411 

412 

413 

414 

415 

416 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 

432 

433 

434 

435 
1 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 
29 


C 

COMMON/ALPHA/A, B,C,D,E,F,R,CC,SB,CD, AT, VU,TCRIT,TR,PR, 

1 OVTRPR.AS 
C 

C THIS ROUTINE USES NEWTONS METHOD TO FIND THE ROOTS OF THE 
C VAPOR- PRES SURE EQUATION; THE SATURATED TEMP CORRESPONDING TO 
C THE GIVEN P 
C 

C AN INITIAL GUESS FOR X ORIGINAL PAGE 19 

C OF POOR QUALITY 

X=560.0 
C 

C USE NEWTONS METHOD UNTIL ERROR IS LESS THAN TSTOL 
C 

DO 75 K= 1,7 

DZDT® (-D(2)*ALOG( 10.0)/(X**2. ) +D(4 ) *ALOG( 10. ) +D(3)/X 
1 -D(5)*D(6)*ALOG(D(6)-X)/(X**2.) +D(5)/X)*EXP(ALOG(10. ) 

1 * (D( 1 )+D(2)/X +D(4)*X) +D(3)*ALOG(X) +D(5)*(D(6)“X)* 

1 ALOG(D(6)-X)/X) 

Z=EXP((D(1)+D(2)/X +D(4)*X)*ALOG(10. ) +D(3)*ALOG(X) + 

1 D(5)*(D(6)-X)*ALOG(D(6)-X)/X) -PRESS 

C 

C COMPUTE NEW VALUE FOR TEMP SATURATED 
C 

X = X— Z /DZDT 

IF(Z/DZDT .LT. ERROR) GO TO 75 
75 CONTINUE 

RETURN 
END 
C 

C i 

C THIS ROUTINE COMPUTES THE MASS FLUX ACROSS THE LIQUID VAOR INTERj 
C FACE, THE EVAPORATION RATE. DUHAMMELS SUPERPOOSITION INTEGRAL IS 
C USED IN APLYING THE SEMI -INFINITE SOLID WITH TRANSIENT SURFACE , 
C TEMPERATURE • 

C IMPROVED MASS USING NEW INDICIES TO GIVE PHI(i)=TS 
C 

SUBROUTI NE MASS ( Y , T , TS , HFGTS , KLTS , VLTSP , MI , ENDRKS , ERF ) 
COMMON/ALPHA/A , B, C,D,E,F,R,CC, SB, CD, AT, UV,TCRIT,TR, PR, 

1 UVTRPR , AS 

REAL ENDRKS 

REAL PHI (6) , KLTS, MI ,ERF( 150,2) ,THETA( 100) 

REAL U(10, 10) ,MT( 10, 10) ,MX( 10) ,MB( 10) ,MR( 10) 

INTEGER N,NN,IV( 10) 

C 

NN* { T+0 . 0 1 ) /0 . 05 + 1 
IF(NN.LT. 2 ) SAVED* 0 . 0 
C 

C NN IS THE NUMBER OF TIME STEPS WHICH HAVE TAKEN PLACE UP TILL NOV* 
C 1 

C 

C COMPUTE MASS FLOW RATE 1 TIME PER RUNGE-KUTTA STEP 
IF ( ENDRKS. EQ. 0.0) GO TO 123 
C 
C 

THETA (NN)«TS 
I F ( NN . LT . 2 ) GO TO 123 
C 



Listing of MAIN* 
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30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 


c 

c 

c 

c 

c 


c 


c 

c 

c 

c 


AS IS THE LIQUID SURFACE AREA IK FT2 
AS=0. 0599332 

CALCULATE SPECIFIC HEAT OF LIQUID BY LINEAR CURVE FIT 
UNITS FT-LBF/SLUG-R 

CLTS=(TS*0. 00003 1666+0. 190 144) *778. 16*32.174 
CALCULATE AALPHA , THERMAL DIFFUSIVITY 
AALPHA=KLTS * VLTSP/CLTS 

COMPUTE THE DEPTH AT WHICH THE TEMPERATURES IN THE FLUID 
WILL BE APPROXIMATED. THE PENETRATION DEPTH IS FOUND, AND THEN - 
OF THIS VALUE IS USED AS THE REG. ON IN WHICH THE TEMPERATURES Wllj 
BE DETERMINED. THIS DEPTH IS THEN DIVIDED INTO 6 LOCATIONS. j 

DEPTH=0. 10*1 .39*2. *( (AALPHA*T)**0.5)/6. 

i 

COMPUTE TEMP AT SIX LOCATIONS, STARTING AT THE LIGUID-VAPOR IN- 
TERFACE USING DUHAMMEL’S SUPERPOSITION APPLIED TO A SEMI-INFINITE 
SOLID WITH TRANSIENT SURFACE TEMPERATURE 
DO 88 I “1,6 
DELX=(I-1)* DEPTH 


IF TIME“0 . 0 , LIQUID IS UNIFORM TEMP AT TSAT 


IF(T.EQ.O.O) PHI ( I ) “TS 
IF(T.EQ.O.O) GO TO 88 

PHI ( I ) «THETA ( 1 ). 

DO 90 K= 2 , NN 
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DELT“T- (K-2 )*0 . 05 


VAL=DELX/ { 2 . * < DELT* AALPHA ) * * 0 . 5 1 
FIND ERF(VAL) 


DO 77 J-1, 102 

I F ( VAL . LT . ERF ( J , 1 ) ) GO TO 83 
77 CONTINUE 

WRITE(6,5) J 
5 FORMAT (13) 

83 ERFVAL=ERF{ J- 1 , 2 ) +(ERF( J , 2 ) -ERF ( J- 1 ,2) )* (VAL~ERF( J-1 , 1 ) )/ ■ 

1 (ERF ( J , 1 ) -ERF ( J- 1,1)) 

19 ERFC“ 1 -ERFVAL 

PHI ( I ) « ( THETA ( K ) -THETA ( K- 1 ) ) *ERFC+PHI ( I ) 

I F ( K . GT . 7 0 ) GO TO 90 
90 CONTINUE 

88 CONTINUE 

SET UP THE COEFFICIENT MATRIX FOR A LEAST SQUARES THIRD 
ORDER CURVE FIT. 


N*4 

U( 1 , 1 )-6. 



L‘i stint) of MAIN+... 
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88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 


0(1,2)«15.*DEPTH 
U(2, 1)=15.*DEPTH 
U( 1 , 3)=55.*DEPTH**2. 

0(2 ,2)=55.*DEPTH**2. 

U(3, 1)=55.*DEPTH**2. 

0(1 ,4)=225.*DEPTH**3. ORIGINAL PAGE IS 

U(2, 3)=225. *DEPTH**3. OF POOR QUALITY 

U(3, 2)=225. *DEPTH*#3. 

U(4,1)=225.*DEPTH**3. 

0(2,4)=979.*DEPTH**4. 

0(3, 3)=979. *DEPTH**4. 

0(4 , 2)=979. *DEPTH**4. 

0(3,4)=4425.*DEPTH**5. 

0(4 , 3)“4425.*DEPTH**5. 

O(4,4)=20515.*DEPTH**6. 

C 

C IN ORDER FOR MORE ACCURATE MARTRIX ARITHMETIC, THE VALDES OF PHI 

C WILL BE SCALED DOWN TO THE SAME ORDER OF MAGNITUDE AS THAT Of 

C DELTA X. 

C 

C WRITE (6, 65) THETA ( 1 ) 

C WRITE(6, 65) TS 

DO 66 11=1,6 

PHI ( I I ) =PHI ( I I ) -THETA ( 1 ) 

C WRITE(6,65) PHI (II) 

65 FORMAT ( F 1 5 . 6 ) 

66 CONTINUE 
C 

C 

MB( 1 )=PHI ( 1 )+PHI (2)+PHI (3)+PHI (4 )+PHI (5)+PHI (6) 

MB ( 2 ) =DEPTH* ( PHI (2 )+2*PHI ( 3 )+3*PHI (4)+4*PHI (5) 

1 +5*PHI (6) ) 

MB{ 3) “DEPTH* *2* (PHI (2)+4*PHI (3) + 9*PHI (4)+16*PHI (5) 

1 +25*PHI (6) ) 

MB (4 ) “DEPTH* *3* (PHI ( 2 ) +8* PHI ( 3 ) +27*PHI (4 ) +64*PHI(5) 

1 +125*PHI (6) ) 

C 

C CALL THE SUBROUTINES SLUD AND SLIR. SLUD COMPUTES THE LO-DECOMP- 

C OSITION OF THE MATRIX U. SLIR COMPUTES A SOLUTION TO THE SYSTEM 

C OF LINEAR EQUATIONS U*MX=MB. 

C 

CALL SLUD(N, 10,U, 10,MT,IV) 

CALL SIR(N, 10, U, 1 0 ,MT, IV,MX,MB,MR, IER) 

C 

C 

C WRITE(6, 122) MX(1),MX(2),MX(3),MX(4) 

122 FORMAT (4E 15. 8) 

C 

C 

DTDX0“MX(2) 

SAVED® DTDXO 
GO TO 124 

123 DTDX0= SAVED 

C ) 

C COMUTE MASS FLOW RATE BASED ON THE SLOPE AT THE INTERFACE, REPRE-j 
C SENTED BY MX(2). 

C 


124 


MI =AS*KLTS* DTDX 0 /HFGTS 
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146 RETURN 

147 END 
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1 

2 

INITIAL VAPOR TEMPERATURE 
INITIAL VAPOR MASS IS 0. 

IS 531. 
23703E-02 

740RANKINE 

LBM 

ORIGINAL PAGE i$ 

3 

NOZZLE AREA IS 

0.31490E- 

04 SQUARE 

FEET 

OF POOR QUALITY 

4 

5 

6 

ULLAGE VOLUME IS 0.68154E-02 CUBIC PEET 

DISCHARGE COEFFICIENT IS 0.77 

TIME T VAPOR TSAT P VAPOR VAPOR MASS 

VENT RATE 

EVAP RAT 

7 

8 
9 

SECS 

KELVIN 

KELVIN 

PASCALS 

KG 

KG/SEC 

KG/SE 

0.05 

295.3511 

294.1243 

91584.75 

0.00102505 

0.9752E-03 

0.0 

10 

0.10 

295.3027 

292.8953 

87597.38 

0.00097898 

0. 9328E-03 

0.3253B-D 

11 

0.15 

295.2554 

291.7068 

83870.31 

0.00093605 

0.8932E-03 

0.5444E-Q 

12 

0.20 

295.2085 

290.5537 

80373.00 

0.00089588 

0.8561E-03 

0. 7133E-Q 

13 

0.25 

295.1621 

289.4314 

77084.81 

0.00065822 

0.821 IE-03 

0.8537E-0 

14 

0.30 

295.1165 

288.3411 

73989.19 

0.00082285 

0 . 7882E-03 

0.9737E-C 

15 

0.35 

295.0713 

287.2834 

71085.00 

0.00078975 

0.7574E-03 

0.1108E-d 

16 

0.40 

295.0264 

286.2517 

68333.63 

0.00075845 

0 . 728 IE-03 

0.1170E-a 

17 

0.45 

294.9817 

285.2500 

65745.81 

0.00072908 

0.7006E-03 

0.1270E-C1 

18 

0.50 

294.9373 

284.2720 

63295.86 

0.00070133 

0.6746E-03 

0.1327E-C 

19 

0.55 

294.8931 

283.3237 

60991.94 

0.00067529 

0. 6501E-03 

0. 1415E-C 

20 

0.60 

294.8489 

282.4033 

58814.86 

0.00065071 

0 . 6269E-03 

0.1472E-G 

21 

0.65 

294.8049 

281.5078 

56759. 19 

0.00062755 

0 . 6051E-03 

0. 1528E-C 

22 

0.70 

294.7610 

280.6326 

54806.49 

0.00060556 

0 . 3843E-03 

0. 1554E-0 

23 

0.75 

294.7170 

279.7825 

52959.86 

0.00058484 

0.5647E-03 

0. 1597E-C 

24 

0.80 

294.6731 

278.9617 

51223.87 

0.00056536 

0 . 5462E-03 

0. 1661E-C 

25 

0.85 

294.6292 

278.1584 

49571.90 

0.00054686 

0.5286E-03 

0. 1674E-q 

26 

0.90 

294.5852 

277.3801 

48010.46 

0.00052938 

0.5120E-03 

0. 1710E-Q 

27 

0.95 

294.5415 

276.6360 

46552.53 

0.00051309 

0 . 4965E-03 

0.17B5E-C 

28 

1.00 

294.4976 

275.9126 

45170.95 

0.00049767 

0.481 8E-03 

0. 180BE-Q 

29 

1.05 

294.4534 

275.2104 

43861.89 

0.00048307 

0.4679E-03 

0.1830E-CI 

30 

1.10 

294.4092 

274.5315 

42624.61 

0.00046929 

0.4548E-03 

0.1B58E-0j 

31 

1. 15 

294.3645 

273.8606 

41428.80 

0.00045598 

0.4420E-03 

0. 1823E-0I 

32 

1.20 

294.3198 

273.2078 

40293.94 

0.00044336 

0.4300E-03 

0. 1837E-0 

33 

1.25 

294.2751 

272.5894 

39239.15 

0.00043164 

0.4187E-03 

0.1 90 1E-0 

34 

1.30 

294.2300 

271.9839 

38228.84 

0.00042042 

0.4080E-03 

0. 1892E-0 

35 

1.35 

294.1846 

271.3838 

37248.34 

0.00040955 

0.3976E-03 

0. 1854E-0) 

36 

1.40 

294.1389 

270.8091 

36328.45 

0.00039935 

0.3878E-03 

0. 1B89E-0 

37 

1.45 

294.0933 

270.2588 

35464.63 

0.00038979 

0 . 37B6E-03 

0. 1920E-0 

38 

1.50 

294.0471 

269.7190 

34633.59 

0.00038059 

0 . 3698E-03 

0. 1903E-0 

39 

1.55 

294.0010 

269.2026 

33853.73 

0.00037196 

0.3615E-03 

0. 1932E-0 

40 

1.60 

293.9543 

268.7019 

33110.49 

0.00036375 

0 . 3536E-03 

0. 1933E-0 

41 

1.65 

293.9075 

266.2014 

32380.82 

0.00035569 

0.3458E-03 

0. 1886E-0 

42 

1.70 

293.8604 

267.7234 

31696.46 

0.00034814 

0.3385E-03 

0. 1912E-0; 

43 

1.75 

293.8127 

267.2688 

31055.70 

0.00034107 

0.3317E-03 

0. 193BE-0 

44 

1.80 

293.7651 

266.8276 

30445.05 

0.00033434 

0 . 3252E-03 

0. 1939E-0' 

45 

1.85 

293.7170 

266.3921 

29850.50 

0.00032779 

0.3189E-03 

0.191 1E-0' 

46 

1.90 

293.6687 

265.9607 

29271.58 

0.00032141 

0.3127E-03 

0. 1BB4E-0 

47 

1.95 

293.6201 

265.5596 

28741.83 

0.00031559 

0 . 307 IE-03 

0. 1934E-0 

48 

2.00 

293.5713 

265. 1802 

28247.04 

0.00031014 

0.301 8E-03 

0. l957E-0i 

4.9 

2.05 

293.5220 

264.7917 

27748.74 

0.00030467 

0 . 2965E-03 

0. 1897E-0j 

50 

2.10 

293.4724 

264.4243 

27283.37 

0.00029955 

0.2916E-03 

0. 1919E-0 

51 

2.15 

293.4226 

264.0605 

26828.18 

0.00029455 

0.2868E-03 

0. 1B92E-0' 

52 

2.20 

293.3726 

263.7139 

26401.30 

0.00028987 

0 . 2622E-03 

0. 1909E-0I 

53 

2.25 

293.3223 

263.3708 

25983.02 

0.00028528 

0 . 2778E-03 

0. 18B3E-0; 

54 

2.30 

293.2715 

263.0532 

25601.21 

0.00028110 

0.2737E-03 

0. 1921E-0 

55 

2.35 

293.2205 

262.7253 

25212.00 

0.00027683 

0 . 2696E-03 

0. 1864E-0 

56 

2.40 

293.1692 

262.4270 

24861.89 

0.00027300 

0.2659E-03 

0. 191 1E-0 

57 

2.45 

293.1177 

262.1365 

24524.09 

0.00026930 

0 . 2623E-03 

0. 1902E-0 

58 

2.50 

293.0659 

261.8228 

24164.07 

0.00026536 

0.2584E-03 

0. 1815E-0 
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59 

2.55 

293.0137 

261.5383 

23841.94 

0.00026184 

0.2550E-03 

0.1863E-C 

60 

2.60 

292.9614 

261.2417 

23509.52 

0.00025820 

0.2515E-03 

0.1806E-C 

61 

2.65 

292.9087 

260.9753 

23214.31 

0.00025497 

0.2484E-03 

0. 1855E-( 

62 

2.70 

292.8560 

260.7061 

22918.13 

0.00025174 

0.2452E-03 

0. 1821E-( 

63 

2.75 

292.8027 

260.4424 

22631.80 

0.00024861 

0.2422E-03 

0. 1812E-( 

64 

2.80 

292.7493 

260.1946 

22365.78 

0.00024571 

0.2394E-03 

0. 1828E-( 

65 

2.85 

292.6956 

259.9565 

22111.23 

0.00024293 

0.2367E-03 

0. 1825E-( 

66 

2.90 

292.6418 

259.7234 

21864.62 

0.00024025 

0.2340E-03 

0. 1816E-( 

67 

2.95 

292.5876 

259.5034 

21635.49 

0.00023775 

0.2316E-03 

0. 1830E-C 

68 

3.00 

292.5334 

259.2656 

21388.68 

0.00023506 

0.2290E-03 

0. 1765E-C 
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- r *~ttXS* uni 
1 

• 2 * 

3 

4 

j 

V 

8 

? 

10 

11 

u 

13 

14 
13 
It 

17 

18 
1 ? 

21 * 

2 i 
* i 
l* 
i n 

IZ 

k w 
2 ? 
it 

J f. 


34 

35 
34 
37 
28 
5 ¥ 

«U 

41 

42 

42 
44 

43 

44 

47 

48 
4¥ 

58 

51 

52 
2 2 
54 
5 3 
54 

57 

58 
5 ¥ 

60 
4 1 
4 2 

43 

44 

• 4 

* w 

44 
4 7 
48 
47 
7lr 

71 

72 
7 2 
74 
VS 
74 

7 7 
78 
77 
8U 

8 i 
(2 
82 

-* 2 _. 


* 


80.00. 

8 82.0. 02254 , 
8 04,0.04511, 
v 84,2.04742, 
0 08,007003, 
2 20.0.11244, 
0 12,0.13474 . 
8 14, 0.15475. 
0 14,0.17701, 
8 18,0.20074, 
8 20.0.22270, 
0.22,0.24430, 
0 24,0.24370. 
8 24.028470, 
0 28,0.30788, 
8 . 38,0 32C63, 

0 32.0.34713. 
8.34,8.34734, 
834.0 30733. 
2 38 .8 40701 . 
U 40,0.42837, 
2 42.8 44747, 

1 44.8 44422, 
v 44,2-48444, 
8 48.0 50275. 

2 58,0.52030. 
1-52.8.33770 . 
2 54,0.55474. 
8 54,8.57142, 
830.058772, 
b 6 8,0. 48384 . 
2 42,8.41741, 
8 44,8. 434S7 , 

8. 44. 0. 44 V 38, 
8 68,8.44278, 
2 .'8,8.47780, 
8 72,8.47143, 
i-.74,U 78448 . 
8 74,0.71754, 

8.78.0. 73001 , 
0.80,0. 74210 , 
8 82,2.75381, 
8 84,0.74514. 

084.0. 77410, 
0.88,0 78449 , 
0 VC. 0.77491, 

0 72,0.80477, 
8 74,8.81427, 
8.74,8.02542, 

8.78.0. 83423, 

1 88,8 84270, 
1 22,8.85084, 
1 84,8 85845. 
1 84.0 86414, 
1 08,0 87333 , 
1.18.8. 88020, 
1 i 2 , 8 88079, 
1 14,8 8 V 3 0 8 , 
1 1 4.8 07710 , 
1 18,0 78484, 
1 20,0 91031. 
1 62,2 VI 53 3, 
1 .24 , 0. 92858 . 
1 26,2.72524, 
1 28,8 727 73 . 
1 38,2.73401, 
1 32,8. 73884 , 
1 34,8 .74171 , 
1 36.8 74554. 
1 28,8.74902, 
1 48,8 ¥3228, 
1 42,8.73538, 
1 44, U 75838, 
. 4* ,8 76103, 
1 48,0 76365, 
1 50.8 76610, 
1 52,0. 94841 , 
1 562 77263, 
1 68,8 9/435, 
1 64,0 77762, 
, 68,8 98247, 
1 72,878580. 
i 74,0 98717, 
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hm*l ... 

v t 

tv 

V9 

' VO • 
VI 
72 
V3 
90 
VS 
Vi 
V7 
va 

V9 

100 

101 

102 


1 96.0.99003, 

2 00, 0 V93322. 

210.0 VV7020 , 

2 . 20.0. 998137, 

2.30.0, 998857, 

2 .00. 0.999311, 

2 .50.0. 999593, 

2 .60.0. 999730, 
2 70,0 V V V 0 6 6 , 

2.80.0 . 999923, . 
2 Vo , 0 . 999939. 

3 .00. 0.9 99978, 

3 .20.0. 999990, 

3. 00. 0.999998, 

360.1 000000 , 
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ALGORITHM FOR DERIVING T> 


1 1 


YNAMIC PROPERnES 


Four governing equations were obtained from DuPont, (Ref. 3) for R-ll; the 
equation of state, the vapor pressure curve, density of the saturated liquid, 
and the heat capacity of the vapor. 

P=P(v,T); equation of state. 

Psat=Psat(Tsat); vapor-pressure curve. 

ft * U C&A* density of saturated liquid. 

Csi'-Cv" (t)' heat capacity of vapor. 

From these four equations, and given Tv and mv, the thermodynamic properties of 
the liquid and vapor may be determined as follows. Refer to Fig. ci, a T-S 
diagram, to Identify the states oelng determined. 

1) • % 

2) P=P(v,Tv); determine system pressure from equation of state. 


3) T$at*T$at(°); determine Tsat from vapor-pressure curve. 


4) Find *tC(Tr,Pr), *i£(Tsat,P), *t£(T$at, Pr), Tv, Pr), from equation of 
state. 


5) Cv «Cv (Tv); find heat capacity of vapoi . 

<i/,(T„.P) 


6) u (Tv, P") = (Zrg'Pji^i f Cv„ iT f t U (Tr Fr) 
Mr*, Pr) T r 




5%l paqe ® 

POOR QUALITY 





A°PENDIX D 


ORIGINAL PAGE f9 
OF POOR QUALfTY 


ADIABATIC DECOMPRESSION MODEL 


In order to evaluate the effect of interfacial mass transfer on the ullage 
pressure response, an adiabatic model was constructed. Derivation is identical 
to that of the interfacial mass transfer model discussed in the ANALYSIS section 
of this report, except that the evaporation rate, ml, is assumed to be zero. 

The continuity equation, Eq. ft), then becomes 


= (Dl) 


The energy equation, Eq.(8), becomes 

(D2) 

These two equations, combined with Eq.(9) now define the vapor space behavior. 
The computer algorithm in App. B is easily modified to solve these governing 
equations. The MASS subroutine, which calculates ml. Is removed and in place is 
put ml =0.0. The remainder of the program is unchanged. 
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APPEhOIX E 

ANALYSIS OF PAST VENTING MODELS 


As discussed earlier, the critical element in modeling the pressure 
response of a cylinder initially filled with a saturated mixture and slowly 
vented is the method used to evaluate interfacial mass transfer. Labus, et al 
(Ref. 1) used the equation 


hht — Ai fC* ( ~ 1 c ) (El 

(poctY'*- h + j 

This equation was obtained by simplifying an analytical expression for the 
interfacial mass transfer during depressurization for an infinitely planar 
interface obtained by Thomas and Morse (Ref. 8). Now, assuming that there is 
no heat transfer across the interface, equations (11) and (12) again apply 



- L A — I 

“ * * t A* U 


d* U *0 


(E2) 


With the definition a=k/pc, equation (El) becomes 


* _ A; A (To m Li ) 

i ■** • — 

Now, solvin; far dT/dxj^ from equation (E2): 


(E3) 


-iC/ _ ( T 0 - Ti ) (E4) 

4K\k-o ( frctV 7 *- 

Equation (E4) is tne temperature gradient of the liquid at the interface, and 
is precisely the temperature gradient at the surface of a semi-infinite planar 
solid undergoing a step change in surface temperature(ref .2) . But, the system 
being modeled undergoes a transient change m surface temperature. Hence, 
some method of incorporating this transient effect, such as Dunammel's 
superposition intearal must be employed for proper applicaio.-. of equation 
(E4). 


Ir, deriving equation (El), labus, et al made a number of assumptions which 
greatly reduced the complexity of the equation derived by* Thomas and Morse. It 
wcs assumed that Tv*Tsat I Pv. The effect of this assumption was discusssed 
earlier. Also, a term in the original expression of Thomas and Morse was 
dropped, assuming the effect of that term to be negligible. The validity of 
this assumption was not evaluated. The equation derived by Thomas and Morse 
wes not used in the present work. Future models may wish to evaluate the 
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